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We study the characterization of multipartite entanglement for the random states of an n-qbit 
system. Unable to solve the problem exactly we generalize it, changing complex numbers into real 
vectors with N c components (the original problem is recovered for N c = 2). Studying the leading 
diagrams in the large-7V c approximation, we unearth the presence of a phase transition and, in an 
explicit example, show that the so-called entanglement frustration disappears in the large- N c limit. 
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I. INTRODUCTION 

The study of entanglement is almost as old as quan¬ 
tum mechanics, as it was the subject of seminal papers 
by Einstein, Podolsky and Rosen [1] and Schrodinger [2]. 
The focus of the founding fathers was on the puzzling, 
non-classical aspects of quantum correlations. Nowadays, 
entanglement [3-5] is viewed mostly as a crucial resource 
in quantum applications, quantum communication and 
quantum information processing [6], and it mantains its 
original fascination for the comprehension of the struc¬ 
ture and geometry of quantum mechanics [7]. 

While bipartite entanglement is well understood and 
quantified, the notion of multipartite entanglement is 
more elusive. This is due to a number of concomitant 
factors. First of all, for many-body quantum systems, 
the number of entanglement measures grows exponen¬ 
tially with the system size [8], making a characterization 
of quantum correlations complicated [9]. Second, new 
properties arise when more quantum parties are involved, 
among these the intriguing appearance of frustration. In 
agreement with the classical notion [10], this is related to 
the impossibility of satisfying a number of requirements 
at the same time [11]. Applied to entanglement, this 
means that given three (or more) parties A, B and C, if 
the entanglement between A and B grows, that between 
A and C or B and C decreases. This is also referred to 
as “monogamy” of entanglement [12- 14]. 

In our work the symptoms of such frustration will sur¬ 
face in the investigation of n = 4 qubits [15, 16] and for 
small number of qudits [17-21], and cannot be avoided 
for n > 8 qubits, where one can prove that it is impossi¬ 
ble to maximize the bipartite entanglement by maximiz¬ 
ing indipenclently all possible bipartitions of a system of 
qubits [22-24]. Further interesting contributions to these 
problems were given in Refs. [25, 26]. 


In this article we will study the properties of multi¬ 
partite entanglement by adopting the concepts and tools 
of classical statistical mechanics [27-30]. This approach 
has proved useful in the study of the bipartite entan¬ 
glement of a large number of qubits, where it unearths 
the presence of phase transitions [31, 32]. The situation 
with multipartite entanglement is, however, much more 
complex. This is due to the fact that the monogamy of 
entanglement acts effectively as a frustration, and the re¬ 
sulting statistical system is frustrated, to say, similarly 
to a spin glass. 

In order to explore the rich landscape that ensues, we 
shall make use of techniques that are based on the analy¬ 
sis of diagrams that naturally arise when one considers a 
liigh-temperature expansion of the distribution function 
of the measure of multipartite entanglement (the poten¬ 
tial of multipartite entanglement) [28]. 

Unfortunately, the evaluation of the contributions of 
different kinds of graph and their resummation is not 
a simple task. One would like to find a strategy that 
can select a procedure to sum some particular families 
of diagrams in order to analyze the features of multipar¬ 
tite entanglement. In this article we will therefore intro¬ 
duce a new parameter, the “colour” index N c (the case of 
qubits is recovered when N c = 2). This appears in a nat¬ 
ural way by generalizing of the structure of the measure 
used to characterize the multipartite entanglement. We 
will consider the limit of large N c and reveal the presence 
of a phase transition. Moreover, we will see an explicit 
example where the frustration of multipartite entangle¬ 
ment disappears if the value of N c is large enough. This 
simplification for large N c is common in statistical mod¬ 
els and field theories and the behavior of the larg e-N c 
approximation usually captures some (but never all) the 
behavior of the finite, or even small-./V c theories. 

This article is organized as follows. In Sec. II we in- 
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troduce notation and familiarize with the statistical me¬ 
chanics approach we will use. In Sec. Ill we evaluate the 
first relevant cumulants of the theory. We restrict our at¬ 
tention to an interesting class of diagrams in Sec. IV: this 
enables us to look at the limit of large number of colors. 
We investigate the behavior of the energy (representing 
entanglement) in Sec. V: this unveils the presence of a 
phase transition, that is studied in Sec. VI. The phase 
transition is further investigated in Sec. VII, where we 
numerically show that no hysteresis appears. We con¬ 
clude in Sec. VIII with a few comments. 


II. MULTIPARTITE ENTANGLEMENT AND 
STATISTICAL MECHANICS 


A. Potential of multipartite entanglement 


Let us consider an ensemble S = {1,2, ...,n} of n 
qubits in the Hilbert space Hs = (C 2 )®". In this article 
we will focus on pure states 


I z) = ^2 z k\ k), Z k G C, 

5Z \ z k\ 2 = 1, 

(2.1) 

fee zj 

fee zj 


where z = (z k ), k = ( ki ) ieS , fc* 

GZ 2 = {0,1}, and 


|&> = 0 | ki)i, \kf)i G C 2 

, ( k j kj ) — S’ij . 

(2.2) 


ies 


Consider a bipartition (A, A) of the system, where A C 
S' is a subset of n A qubits and A = S\A its complement, 
with ua + tia = n. We set ua < ua without loss of 
generality. The total Hilbert space factorizes into Hs = 
Ha®Ha, with Ha = Cf, of dimensions Na = 2 nA 

and Na = 2 nA , respectively (Na^A = AT). 

The bipartite entanglement between the two subsets 
can be measured by the purity of the reduced state of 
subsystem A 

n A {z) = Tr p\, p A = Trql'i/’XVi (2.3) 

Tr^ being the partial trace over A. We notice that tta = 
7r a and 


1/N a <n A <l, (2.4) 


where ua = [n/2] (balanced bipartitions), with [a:] being 
the integer part of x. The coupling function is [24, 29] 

A = g{(k®l)V(k'®l'),(k®l')\/{k'®l)), (2.6) 

where 


g(a, b) = 5 aAb , 0 g(\a\, |6|) (2.7) 


and 


g(s,t) 


1 ( n\ fn — s — t\ 

2 \ n A) [\ n A - s J 



In Eqs. (2.6)-(2.7) we have defined |a| = Ylies ai ’ N = 
a © b = (a* + b z mod 2) i6 g being the XOR 
operation, a V b = (a.j + 6,; — aibi)i^s the OR operation, 
a A b = (a,ibi)i e s the AND operation. Due to its linear 
structure, H(z) inherits the upper and lower bound of 
the purity n A (z) in (2.4): 


1/N A < H{z) < 1. (2.9) 


We notice the following symmetries of the coupling 
function: 


[ = A (k',h,l,l') 

\ A(fc, k'\l, l') = A(Z, l'-, k, k') , (2.10) 

( A(k,k'-,l,l') = A (k',k;l',l) 

which reflect the reality of H(z) 


H{z) = H(z) = H{z). (2.11) 


Moreover, since 

A(fc © to, k 1 © to; l © to, l' © to) = A(fc, k 1 ; l, l'), (2.12) 

the Hamiltonian is invariant under rotations and reflec¬ 
tions: 


H(z mm ) = H(z k ). (2.13) 


for every to G Iff. 


B. Classical statistical mechanics approach 


where the lower and upper bounds are obtained for max¬ 
imally entangled and separable states, respectively. 

A natural extension of this measure to the multipartite 
scenario is realized considering the potential of multipar¬ 
tite entanglement , that is the average bipartite entangle¬ 
ment between balanced bipartitions [29]: 


The analysis of the properties of H(z) can be rephrased 
in a classical statistical mechanical framework. Let us 
consider the partition function of a system with Hamil¬ 
tonian H{z) at a fictitious temperature /3 _1 , 

Z(p) = J e-tnW d p(z), (2.14) 

where [33] 


= (nj £ ^(*) 

V 7 |A|=n A 

= ^2 A (k,k';l,l')z k Zk'ZiZi', (2.5) 

k,k',l,l'e zj 


d n(z) = N ^' <5(1 - ||z|| 2 ) \\&z k &.z k , 

k 


(2.15) 
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is the uniform measure on the hypersphere 

ni 3 = Ew 2 = 1 . ( 2 - 16 ) 

k 

and dzkdzk denotes the Lebesgue measure on C. 

The Lagrange multiplier p fixes the average value of 
multipartite entanglement. In particular, [3 = 0 corre¬ 
sponds to the uniform sampling of random states. The 
limits p —> ±oo select configurations that, respectively, 
minimize or maximize H(z)\ the former case corresponds 
to maximally multipartite entangled states (MMES) [23], 
the latter to completely separable states. 

The average value of H (z) (seen as the average energy 
of the system) at arbitrary ft can be obtained as 

{H)p= zW)J He ~ PH d^=-^ lnZ ^)- 

The properties of the distribution function of the po¬ 
tential of multipartite entanglement can be analyzed by 
evaluating its cumulants. In particular, the m-th cumu- 
lant of H(z) reads 

f)m r\m— 1 

4 = 

(2.18) 

III. DIAGRAMMATIC EVALUATION OF THE 
CUMULANTS 



FIG. 1: (Color online) Graph contributing to the first cumu- 
lant: two-leaf cactus diagram. 


numbers defining the wave functions to real vectors of 
N c components and take N c large. In the following two 
sections we recall the computations in [29], and we an¬ 
ticipate which diagrams are going to be relevant in the 
large N c limit. 


A. First cumulant 

The average energy at p = 0 reads 

(H) 0 = E ^■( k i^ k 2'ih,l2)( z kiZk 2 Zl 1 Zl 2 )o 

= <ki| 2 |^ 2 | 2 )o E bd)iP( 2 )] (3-3) 

p65 2 

where 


In this section we will briefly review some properties of 
the high temperature expansion of the distribution func¬ 
tion of the potential of multipartite entanglement. We 
remind that for P = 0 one gets the unbiased random 
states. The average energy reads 

/m _ y' to I', 

{H)p ^ (m, — 1 )! K ° (3 ' 1} 

where the brackets (• • • )o denote the average with re¬ 
spect to the uniform (unitarily invariant) measure (2.15). 
The only nonvanishing correlation functions are of the 
form [29] 


[p(l)p(2)}= E A(fci,fc 2 ; fcp(i),fcp( 2 )) (3.4) 

fcl ,£ 26^2 


and S 2 is the symmetric group of order 2. This is repre¬ 
sented by the doubly degenerate graph in Fig. 1 which, 
by (2.6) and (3.2) gives 


(H) 0 = 


Na + N A 
N+l 


(3.5) 


For balanced bipartitions of an even number of qubits, 
Na = Na = y/N, and N -A +00 we get 



(N- 1 + EJLi^)! ’ 


(3.2) 


where rrij are nonnegative integers. In order to calculate 
the required cumulants, we will make use of the diagram¬ 
matic technique developed in [29]. In this way, we will 
be able to evaluate the contribution arising from differ¬ 
ent kinds of graphs. The objective of this study is to try 
to understand whether by embedding this problem in a 
larger family, dependent on a parameter to be introduced 
in the following, called N c , the number of colours, a class 
of graphs can be isolated that dominate in an appropri¬ 
ate limit. To this extent, we will promote the complex 


(H) 0 


2 y/N 
N + l 


2 

Vn' 


(3.6) 


For an odd number of qubits the value of (H )0 is 3/y/2 
larger than (3.6). 


B. Second cumulant 

The second cumulant is defined as 

«o 2) = (H 2 ) 0 - (B)l 


(3.7) 
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FIG. 2: (Color online) Connected graph contributing to the 

second cumulant: three-leaf cactus diagram. FIG. 3: (Color online) Connected graph without leaves con¬ 

tributing to the second cumulant. 



We have 

(H 2 )o = A(fci, h, ^)A(fc3, k±\I 3 , li) 

k,l£l.* n 

x (z kl z k2 z k3 z ki Z h Zi 2 Zi 3 Z U ) 0 

= (N 2 M 2 M 2 M 2 )o [p(l)p(2)^(3)^(4)]• 

p£S 4 

(3.8) 

Let us evaluate the contribution from connected 
graphs. The graph with two links between left and right 
pairs in Fig. 2 has degeneracy 16 and reads 


[1 3,2 4] 


N ( N a + N A f 
4 


(3.9) 


The graph shown in Fig. 3 has degeneracy 4. The asso¬ 
ciated contribution does not have a transparent form; on 
the other hand, its asymptotic formula reads [29] 


and we have 

(H 3 )o = (l-^i | 2 | 2 1 - 231 2 1 ^ 41 “ I ^5 1 2 !^ | 2 )o 

x t^( 1 ) iP( 2 )^( 3 ) ^( 4 )^(5) ^>( 6 )]- 

p£S 6 

(3.14) 

The contributions from connected graphs with three 
leaves (degeneracy 128), represented in Fig. 4(a), 
and with two leaves (degeneracy 192), represented in 
Fig. 4(b), are equal 


[1 6 ,3 2,5 4] = [1 3,2 5,4 6 ] = N ^ Na + Na ^ ~ 7V 5/2 , 

(3.15) 

The graph represented in Fig. 4(c) (degeneracy 192) gives 
an asymptotic contribution 


[3 4,12]- N a (3.10) 

with a = log 2 3 ~ 1.5850. 

Notice the presence of the irrational exponent in the 
graph (3.10). Moreover, since (3.9) is exactly canceled by 
the non-connected contribution from the square of (3.5), 

the graph (3.10) represents the dominant contribution to 
( 2 ) 

the second cumulant k 0 ; that therefore has the asymp¬ 
totic value of 


_(2) _ 3a/2 

K 0 — jy 2.4150... • 


(3.11) 


However, considering only the contribution from the 
graph (3.9), which will turn out to be dominant in the 
large N c limit (see later in the paper), we obtain 


[1 6,2 5,3 4] - 3V2N a+1/2 . (3.16) 

Finally, the asymptotic contributions from the graphs in 
Fig.5(a) and (b) read, respectively [29] 

[5 6,1 2,3 4] ~ N a , (3.17) 

[3 6,5 2,1 4] ~ ciV 7 , (3.18) 

with c ~ 1.05385 and 7 ~ 1.8417. Notice the appearance 
of a new irrational exponent 7 , which again, by cancella¬ 
tions, turns out to be the dominant contribution to the 
cumulant Kq , followed by that due to a. In summary, 
the asymptotic value is 

4 3) ^ 67.4 IV- 4 - 158 -. (3.19) 


r. (2 ) _ 
A, n — 


4 (N a + N A f 2 s 

(N + 1)(N + 2)(N + 3) ~ N 2 


(3.12) 


where the asymptotic expression is valid for balanced bi¬ 
partitions of an even number of qubits, Na = Na = y/N, 
and N —► + 00 . 


C. Third cumulant 


Again, by retaining only the contributions from 
Eq. (3.15), as for Rq , which dominate the large N c limit 
we find 


-.(3) 


40(A)4 + N a ) 


^3 


(N+l)(N- 

320 


2)(N + 3){N- 


NV 2 


•4)(JV + 5) 

(3.20) 


The third cumulant reads 

4 3) = o) 3 ) q 

= ( H 3 ) 0 ~3(H 2 ) 0 (H) 0 + 2(H)l (3.13) 


where the asymptotic expression is obtained for balanced 
bipartitions of an even number of qubits, Na = Na = 
y/Tf, and N —> + 00 . 

For the fourth and fifth cumulant we have to resort 
to numerical methods. By generating a large number of 
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FIG. 4: (Color online) Connected graphs (with leaves) con¬ 
tributing to the second cumulant. (a): four-leaf cactus dia¬ 
gram. (b): four-leaf cactus diagram, (c): one leaf 



FIG. 6: The scaling exponent of K n from numerics. On the 
first 3 points the error bars are smaller than the symbol size 
and the data agrees with the theory. 



FIG. 5: (Color online) Connected graphs (without leaves) 
contributing to the second cumulant. (a): same internal and 
external orientation of the edges, (b): opposite internal and 
external orientation of the edges. 


random vectors for different N, we find that 



= (43 ± 12) at- 418±0 06 

(3.21) 

ft (4) 

= (27 ± 20)IV -5 ' 2±o ' 5 

(3.22) 

ft (5) 

= (148 ± 90)V _6 ' 5±1 ' 5 . 

(3.23) 


(3) 

where Kq ' is in good agreement with the theoretical value 
(compare (3.21) with (3.19)) . 

These numbers suggest a scaling of the form Kq ^ ~ 
jy--o.25-i.4n ( w i lose quality can be seen in figure), in par¬ 
tial agreement with the large- N c result is ~ TV 1-1 ' 5 ". 


temperature expansion Eq. (3.12) and Eq. (3.20), hoping 
that this can give some information about the characteri¬ 
zation of multipartite entanglement. However, the choice 
of the cactus (subdominant) diagrams was quite arbitrar¬ 
ily and motivated only by the fact that they yield integer 
exponents. 

Of course, one should choose the family of graphs with 
an objective criterion, not on the basis of ease of compu¬ 
tation. An objective criterion which is common practice 
in field theory is that of generalizing the model to a larger 
symmetry group, from 50(2) to SO(N c ), and then tak¬ 
ing the limit of large number of colours N c [34-36]. This 
will select a family of diagrams, and enable us to com¬ 
pute their amplitudes in a closed form to leading order 
in 1/N C . We will see that these diagrams are exactly the 
cactuses. Moreover, we will see, at the end of the compu¬ 
tation, that a phase transition is obtained at a temper¬ 
ature below which a combination of the rotation group 
SO(N c ) and the symmetric group S n is spontaneously 
broken. 

In order to perform this calculation, we will rewrite 
the expression of the potential of multipartite entangle¬ 
ment in a different form. The Fourier coefficients Zk of a 
quantum state can be written as 


IV. CACTUS DIAGRAMS AND LARGE 
NUMBER OF COLOURS 

**=$* + *$*, (4.1) 

As the reader should have deduced from the previous 
sections, the proliferation of diagrams and the variety of 
the exponents on N appearing prevent one from writing 

a closed form for the correlation functions (ziZj)p. One with /i = 1,2, real numbers. We will consider this 
can then decide to sum a particular family of diagrams, object as a two components vector. In terms of these real 
as we did for the second and third cumulant in the high quantities, we can rewrite the Hamiltonian in Eq. (2.5) 
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in the following form [using the symmetries in Eq. (2.10)] Now, let 


H = 

k,k',l,l'eZ™ 

2 2 
x 2 £ - 53 <!>£${*, 


H,is =1 




/x,z/=1 
N c 




k,k',l,l'ez™ 


H,is=l 


(4.2) 


for N c = 2, and 

A(fc, jfe' ; l , Z') = 2A(fc, fc'; Z, l') - A(Zs, Z; As 7 , Z'). (4.3) 

This expression can be put in a more compact form by 
writing 

AT _ _ 

£ A(fc, fc'; Z, Z 7 )($fc • 4>/)(4> fc / • <I> ; /), (4.4) 

k,k' ,1,1' ^Z 1 ^ 

where, in general, $ fc = (4>J.,..., $^ c ), and the dot de¬ 
notes the scalar product. The index /i, ranging form 1 to 
7V C , will play in the following the role of a color index, 
and we will be interested in the limit N c —> 00 . The nor¬ 
malization of the complex vector z k becomes, for generic 
N c the constraint 

N c 

EE(« 2 - £**•** = !■ (4-5) 

fcezj ^=1 fcgzj 

We will now show that the Hamiltonian (4.4) is posi¬ 
tive for any color N c , and will find a lower bound which 
generalizes the bound on the potential of multipartite en¬ 
tanglement, valid for N c = 2. From (2.5) and the explicit 
expression of the purity across the bipartition (A, A), 

ka(z) = £ 5k A xJk' A ,iJk A ,i A h' A x A z kZ k 'Ziz v , ( 4.6) 

k,k',l,l' 

it is not difficult to express the Hamiltonian as an average 
over bipartitions 


H(*) = 


riA 


-1 


£ h a ($), 


(4.7) 


\A\—n A 


where 

Ha@) = Y 


Y,Y,( 2 ^ A *& A k A n's*Lv A 


k,k' 


-K A kA kA K' A k' A n Ak ' A ), 

(4.8) 


with kj 1 = (ki) i& A- 


Ak A fi,i A v /* . 4 } 


L 1 (f)V 
k A k A k A l A i 


(4.9) 


k A 


and note that A is a symmetric A 4 N c x N A N C matrix, 
namely X T = A. We get 


H a = 


Nr 


EE( 

k A ,l A 


2X 


k A n,l A v 


Nk A fj,,l A vNk A L',l A fiJ , 

(4.10) 


that is 


N r 


H a = ^ [2tr(A T A) - tr(A T A)] , 


(4.11) 


with Y kAflt i A u = X kA ^i AfJ , a symmetric matrix. By 
the Cauchy-Schwarz inequality for the Hilbert-Schmidt 
scalar product we get 

tr(X T y) < tr(A T A ) 1/2 tr(F T y) 1/2 . (4.12) 

But it is easy to see that tr(F T y) = tr(A T A), so that 


H A ($)>^tr(X 2 )>0. 


(4.13) 


By making use of the constraint (4.5) we can estimate 
the positive lower bound as follows. First notice that the 
positive matrix X > 0 has unit trace. Indeed, 

E = c 

k A ,n k A ,k A ,n 

(4.14) 

by (4.5). Therefore, 


tr(A 2 ) > 


1 1 

> 


rankA N C N A 


(4.15) 


But, from the definition (4.9), we get that rankA < N A , 
so that 


H a ($) > 


N c 
2 N a 


(4.16) 


By a straightforward computation, one can check that 
the minimum is attained at 


k VNA 5kA ’ kA ' 


£(<n 2 = 1, 


(4.17) 


where the Kronecker delta is meant to be 1 when k A = 
(kA, 0,..., 0). Notice that (4.17) is the direct gen¬ 
eralization of the maximally bipartite entangled state 
across the bipartition (A, A), with coefficients z k = 

e iak S kA ,k A /VNA- 

By plugging (4.16) into (4.7) we finally get the desired 
lower bound 


min H > 


N c 

2N a ' 


(4.18) 
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k k k k k k k k 


FIG. 7: (Color online) Dyson equation for the propagator in 
the large N c limit. 


the non-zero propagator, obtained by setting /.t = v, k = 
l, as 


or k = (KKb = gT* 

-0A, E A ( fc ’^ 0 GuGkl 

v,l 

(4.25) 


which may not be attained due to frustration among the 
bipartitions. See the discussion in Sec. V. 

However, an interesting simplification occurs in the 
limit N c —> oo. We now introduce the constraint (4.5) by 
means of a Lagrange multiplier A and rescale (for future 
purposes) the inverse fictitious temperature /3 = /3//3q, 
with 


A) 


2 N 2 

N a + N a - 1 


TV 3 / 2 . 


(4.19) 


We are finally left with the modified Hamiltonian 


0H(\) 


k,k',l,l' 



(4.20) 


From Eq. (4.24) we get 

GZ = - J(| E A ( fc ’(4.26) 

I/,/ 

Let us define the quantity 


n c 


Gi = E G\ 


VV 

ll 5 


I/=l 


that, inserted into Eq. (4.26), gives 

gw _ _ \ _ \ _ 

kk N * N \ + \ 

We now sum over the colour index /i and obtain 

G k = 1 1 


N X + $fr £,A(M;M)G« 


(4.27) 


(4.28) 


(4.29) 


The partition function is an integral over <1> and over 
the Lagrange multiplier imposing the constraint so that 


Z = J dA d$ e“^ (A ). 


(4.21) 


The evaluation of Z can be done by expanding it for small 
/? since the A part is quadratic and, after resummation 
of the diagrams, calculating the saddle point in A. 

The saddle point equation in A is 

= 0 (4.22) 

and it is equivalent to the request that 

E<** • = !. ( 4 - 23 ) 

k 

where the average is evaluated using the full partition 
function Z. 

We should find the value of A(/3) that satisfies the con¬ 
straint. For example, for /3 = 0 we have 

G { T V = <W>o = (4.24) 

where /x, v are colour indices; in order to satisfy the con¬ 
straint we need A = 1. 

In the limit N c —» oo the diagrams giving the dominant 
contribution are the cactuses. The solution for 0 > 0 can 
be obtained from the Dyson equation by considering only 


We also notice that the constraint (4.23) implies 

E °k = 1- (4.30) 

fcGZ 2 

For sufficiently small /?, we look for a solution with 
unbroken permutation symmetry (iSjv-symmetric) 

G k = Vfc G ZJ, (4.31) 

so that from Eq. (4.29) we have 

A +!^E A ( fc ’W) = 1- (4-32) 

l 

From Eq. (2.6) and Eq. (4.3) we get 

Ea(M;M) = J2Mo,k®i)- g {k®i,k®i)] 

l l 

= Et 2 9(0,1)-9(1,1)]. (4.33) 

l 

For balanced bipartitions we obtain 

T,9M= NA ^ NA , E^.0 = 1, ( 4 -34) 

l l 

so that 

Ea(M;M) = N A + N A -1 (4-35) 

l 



independent of k, and Eq. (4.32) reads 

A = 1 -ft. (4.36) 

By using this result we get the average purity (again 
retaining only the leading order in 1/N C ) 

( h ) 0 = ^ £ A(m,M)GM" 

N c (N A +N A ~l) Nc_ 

2 N ~ y/N' [ ’ 

independent of [3 which, for N c = 2, gives the correct 
result for the value of the first cumulant in Eq. (3.6) only 
at (3 = 0. This is due to the fact that, to the lowest 
order in /3, sub-leading diagrams in N c are sub-leading 
in N as well. Notice how, in this approximation, the de¬ 
pendence on the temperature has disappeared (a similar 
phenomenon occurs in a matrix model related to spin 
glasses [37]). 

However, f3 = 1 is a critical temperature, as one can 
see that the $ fluctuations become massless. In fact, the 
Lagrange multiplier A is the coefficient of the quadratic 
part of the hamiltonian (4.20) 


d 2 U 

d<4 


<£=0 


1 


N c N6 a0 6 kl A. 


(4.38) 


where [.] denotes the integer part. We encounter the so- 
called frustration of multipartite entanglement. In other 
words the requirement that the bipartite entanglement 
is maximal (minimal purity) for all bipartitions can en¬ 
gender conflicts. Incidentally, we notice that this phe¬ 
nomenon can be found not only for spin systems but also 
for infinite-dimensional systems as in the case of Gaus¬ 
sian states [38]. 



FIG. 8: Dependence of the rescaled minimum Eq TL ' > of the 
Hamiltonian (4.4) on the value of N c for n = 4. 

In order to study this phenomenon, let us define the 
quantity 


Therefore for f3 > 1, A < 0 and for N — > oo we should ex¬ 
pect spontaneous symmetry breaking such that for some 
k ($fc) > 0. The Sn symmetry gets spontaneously bro¬ 
ken. 

This value of (3 C is not in evident agreement with the 
numerics, although a tendency for large N c of developing 
a kink at (3 ~ 1 is noted in the data for n = 3,4,5 and 
7. In particular in Figs. 9, 10, 11 and 12 one notices that 
for small (3 by increasing N c the data move towards the 
large-iV c , small /3, /I-independent result (4.37), while for 
large (3 the flow is reversed, probably asymptoting to the 
large-_/V c minimum value of the Hamiltonian. 

A complete solution of the partition function of the 
quartic Hamiltonian (4.20) and the 1/N C corrections will 
be the subject of future work. 


V. MINIMUM OF THE HAMILTONIAN IN THE 
LARGE-W LIMIT 

In this section we will give an explicit example of what 
happens to the potential of multipartite entanglement for 
increasing values of N c . 

Let us start by considering a collection of qubits. We 
recall that this case corresponds to N c = 2. In particular, 
it is well known [11, 22, 23] that for n = 4 or n > 8, 
the ideal minimum of the Hamiltonian (4.4) cannot be 
reached i.e. 

£<"> = minH> jL = ^ 2 -[f], (5.1) 


rr( n ) _ 

_ Nr 


(n) 


(5.2) 


which represents the rescaled minimum of the Hamilto¬ 
nian (4.4) for different values of N c . In order to under¬ 
stand the consequence of the large N c limit, we have per¬ 
formed a numerical minimization of Eq. (4.4) for n = 4 
and 2 < N c < 8. In Fig. 8 we plot E^ as a func¬ 
tion of N c . It is manifest that by increasing the value of 
the colour parameter frustration disappears. Indeed, for 
N c = 2 we find Eg”'* = 1/3, in agreement with previous 
results [15, 16, 18]. If N c > 4 we have Eq 1 ' 1 = 1/4. Ap¬ 
parently, the ideal minimum is obtained by minimizing 
each term (4.8) in the Hamiltonian separately as if they 
were independent. Therefore for N c > 4 the system is 
unfrustrated. On the other hand, this means that one of 
the most characteristic trait of multipartite entanglement 
cannot be analyzed in the large-A/, limit. 

Incidentally, it is interesting to notice that for n = 4 
the values of the minimum energy follow the law 


E (4) - 
— 


Nr. 


6Ab 


(5.3) 


and that for N c > 4 this expression becomes < 1/4. 


VI. SEARCH FOR THE PHASE TRANSITION 

In Figures 9, 10, 11 and 12 [respectively for n = 3 
(N = 8), n = 4 {N = 16), n = 5 {N = 32), and n = 7 
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FIG. 9: (Color online) Numerical results for n = 3 (IV = 8) 
and N c ranging from 2 to 20. Bo and Hnc are defined in 
Eq. (4.19) and Eq. (6.1). 


FIG. 12: (Color online) Numerical results for n = 7 (N = 128) 
and N c ranging from 2 to 10. Bo and Hnc are defined in 
Eq. (4.19) and Eq. (6.1). 



FIG. 10: (Color online) Numerical results for n = 4 (N = 16) 
and N c ranging from 2 to 20. Bo and Hnc are defined in 
Eq. (4.19) and Eq. (6.1). 



FIG. 11: (Color online) Numerical results for n = 5 (N = 32) 
and N c ranging from 2 to 20. Bo and Hat c are defined in 
Eq. (4.19) and Eq. (6.1). 


(N = 128)] we plot the expectation value (H)p rescaled 
by 


Hnc 


N c (N a + N a - 1) 
2N 


( 6 . 1 ) 


versus P = P/Po for different values of the color num¬ 
ber N c obtained with MonteCarlo simulations. We recall 
that H is defined in Eq. (4.4) and P 0 in Eq. (4.19). 

In all cases we notice that for large values of N c there 
is an inflection point. Moreover, for increasing values of 
N c , the curves tend to become flatter for small values of 
P/Po- This behavior is more evident for n = 4 (N = 16) 
and n = 7 ( N = 128). Finally, we notice that for n = 4 
with N c > 3 and n = 7 with N c = 10 there is an evident 
change in the behavior of the curves for larger values of 
P/Po, apparently absent in the case of n = 5. This could 
be due to the removal of the multipartite entanglement 
frustration in the case n = 4 as detailed in the previous 
section (the frustration phenomenon is not present for 
n = 5, and the case n = 7 is still open). 

The flattening in the region of small 8 = P/Po might be 
interpreted as a confirmation of the presence of a phase 
transition around p = 1 in the limit of large N c . For 
n — 4 (where this effect is more evident) we notice that 
the inflection point of the curves moves from values of 
P ~ 2 to P ~ 1.5. 


VII. SEARCH FOR HYSTERESIS 

We now try and elucidate the features of the phase 
transition brought to light in the preceding section. In 
particular, we will search for the presence of hysteresis 
when the system is “cooled” fixing a large value of p and, 
after equilibration, is taken back to the initial value of p. 
In order to perform this analysis we will use simulated- 
annealing-like algorithm [39, 40]. After fixing an initial 
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20 temperatures 500 MCsteps 




FIG. 13: (Color online) Search for the hysteresis phenomenon. 
Numerical results for n = 4 (N = 16), N c = 20. /?o and 
Hnc are defined in Eq. (4.19) and Eq. (6.1). Different curves 
correspond to different number of MonteCarlo steps between 
the decrease of f3 in the simulated annealing procedure. See 
the text for details. 


value of p we let the system reach the equilibrium. At 
this point we start to “heat” the system and decrease the 
value of P fixing the number of MonteCarlo steps before 
the following decrease. 

In Fig. 13 we show the results obtained for the case 
n = 4, N c = 20. We start from p = 130 (point not 
shown in the figure), corresponding to P/Po — 1.777, 
and anneal for 20 values of P (each step being equal to 
A p = 4). The different curves correspond to 50, 100, 
200 and 300 Monte Carlo steps before each decrease of 
p. Notice that the second procedure is what one calls 
a quench in the literature (at variance with annealing). 
We cannot see any evidence of hysteresis. As a reference, 
we have included in Fig. 13 also the curve (blue dots) 
that we have already included in Fig. 10: we remember 
that each point has been obtained not with simulated 
annealing but starting from p = 0, fixing the value of 
P, let the system reach the equilibration and then per¬ 
forming a Monte Carlo run; the other curves in Fig. 13 
coincide within errors. The annealing was repeated for 
10 and 5 values of p, with no difference and no evidence 
of hysteresis. 

In order to compare the results using different proce¬ 
dures, we have tried to perform a simulated annealing 
procedure also going from p = 0 to larger values of the 
inverse temperature. We have changed the speed of the 
annealing procedure by varying the number of Monte 
Carlo steps between successive temperatures. The sys¬ 
tem is cooled starting from a given temperature, without 
going back to the initial state at every step. We expect 
to find the same results as in Fig. 13. We start from 
P/Po = 0 and then proceed with simulated annealing 
at steps ~ 0.1 up to P/Po = 3 and then at steps ~ 1 
up to P/Po = 10. The rescaled energy for 500 Monte 
Carlo steps between different temperatures is plotted in 


FIG. 14: (Color online) Numerical results for n = 4 (N = 16), 
N c = 20, and 500 MonteCarlo steps, do and Hnc are defined 
in Eq. (4.19) and Eq. (6.1). See the text for details. 


Fig. 14, where we show (blue dots) the curve obtained 
without the simulated annealing and two different simu¬ 
lations with the simulated annealing (yellow and purple 
squares). The results corroborate our previous finding 
and extend the graph to higher values of /? (below the 
frustrated minimum). No difference is observed for 100 
Monte Carlo steps (not shown). We conclude that no 
hysteresis is present. 

As a final test, we have considered the behavior of the 
overlap between configurations of two different Monte 
Carlo simulation run in parallel (labelled 1 and 2, re¬ 
spectively). This quantity is defined as 


( 9 2 ) 


P 




(7.1) 


where the average is performed at a fixed value of p, by 
considering a configuration every ten Monte Carlo steps. 
We have considered the case of 500 Monte Carlo steps 
between different values of p. The results are shown in 
Fig. 15, where each point is obtained averaging over 50 
configurations at fixed p. The vertical axis is rescaled 
by the quantity 1/NN C which corresponds to the value 
of two random independent configurations. There is a 
decrease of the overlap around P/Po — 2.5, which most 
likely corresponds to different replicas freezing in differ¬ 
ent minima. However, these results cannot be considered 
decisive for the presence of a phase transition. 


VIII. CONCLUSIONS 

We have investigated the behavior of the potential of 
multipartite entanglement by generalizing the problem 
to N c colors and taking a large- N c limit. In the analyti¬ 
cal treatment we find that an instability occurs for suffi¬ 
ciently small temperatures that breaks the permutational 
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500 MCsteps 

N N c (q 2 )p 



FIG. 15: (Color online) Overlap of configurations obtained 
from two different Monte Carlo simulations. The number of 
steps between different successive values of /3 is 500. See the 
text for details. 


symmetry between spin configurations. This instability 
generates, in the thermodynamic limit, a second order 
phase transition. 

We have performed numerics for small system sizes 
(n < 7,N C < 20) and observed the signatures of such 
phase transition (Sec. VII). In order to exclude a first 
order transition, we have investigated the hysteresis phe¬ 
nomenon, not finding any signature of it. The quanti¬ 
tative comparison of the analytic, large N c results, with 
the numerics is plagued by large N c corrections (larger, 
the larger is n), but it is to a large extent satisfactory. 

The frustration generated by the entanglement 
monogamy is to some extent similar to what happens 
in frustrated spin systems which have (for low tempera¬ 
tures) a glassy phase. However, it is worth remembering 
that our model does not contain any quenched disorder. 
Taking this into account, a more fit analogy is with con¬ 
figurational glasses [41-43]. 

Our result is of considerable significance for the physics 


of entanglement. Decreasing the temperature, the av¬ 
erage of the potential of multipartite entanglement de¬ 
creases towards its minimum. A smaller value of the po¬ 
tential is representative of a more entangled state. The 
typical states at low temperature are sampled from a 
measure which is not S n symmetric, so the states which 
minimize the potential are necessarily concentrated on 
fewer spin configurations. This is at first thought coun¬ 
terintuitive. One might think that in order to get more 
multipartite entanglement one should mix in more con¬ 
figurations and create a more uniform wave vector on 
the binary hypercube ITf. We show that the opposite 
is true. The states which maximize the multipartite en¬ 
tanglement have a certain degree of breaking of the S n 
symmetry. 

Finally, there are other features of entanglement frus¬ 
tration that have not been analyzed in this article, that 
pertain to the ground states of some quantum many-body 
Hamiltonians [44-49]. As entanglement is a resource for 
quantum computation, it would be desirable to see if one 
can find a quantum Hamiltonian whose ground state is 
highly entangled. In this direction, recent developments 
have shown that some quantum many-body Hamiltoni¬ 
ans [50-52] made from projectors [53, 54] do violate the 
common area law for entanglement in favor of volume 
or almost-volume law. These ground states are likely 
candidates for quantum certificates of difficult quantum 
computation problems. 
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